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Abstract. In finite-size population models, one can derive Fokkcr-Planck equations 

to describe the fluctuations of the species numbers about the deterministic behaviour. 
1} In the steady state of populations comprising two or more species, it is permissible for 

a probability current to flow. In such a case, the system does not relax to equilibrium 
! K but instead reaches a non-equilibrium steady state. In a two-species model, these 

currents form cycles (e.g., ellipses) in probability space. We investigate the conditions 
& under which such currents are solely responsible for macroscopically-observable cycles 

in population abundances. We find that this can be achieved when the deterministic 

limit yields a circular neutrally-stable manifold. We further discuss the efficacy of one- 
i dimensional approximation to the diffusion on the manifold, and obtain estimates for 

the macroscopically-observable current around this manifold by appealing to Kramers' 

escape rate theory. 
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1. Introduction 

When deciding how to construct an ecological model to describe the dynamics of a 
population, a parameter that controls the population size, called the carrying capacity, 
K, must be taken into account pQ. The traditional approach is to assume that K 
is sufficiently large that fluctuations can be neglected. Then, a population level model 
(PLM) describing the dynamics with ordinary differential equations suffices [2j. However 
for many biological and ecological systems where populations number typically only 
hundreds (as opposed to millions, let alone the vast numbers of particles present in 
condensed matter systems), it is more appropriate to employ an individual based model 
(IBM). Here, the dynamics are cast as stochastic birth and death processes which can be 
described as a continuous time Markov process in a discrete configuration space. This 
configuration space is defined by vector a n = {m, 112, ■ ■ ■ , n^} where n, is the number 
of each of the N species in the population, or equivalently by x = {xi, x 2 , ■ ■ ■ , x^} where 
xi = ni/K. 

The starting point in the analysis of an IBM is the master equation 
dP(x) 



dt 



^[T(x|x')P(x')-T(x'|x)P(x)] (I) 



describing the time evolution of the probability P(x, t) in the space of configurations. 
The birth and death processes are expressed in terms of the transition rates T(x'|x) 
from configuration x to x'. By carrying out a systematic system size expansion in K, 
first prescribed by van Kampen [3], of the master equation the deterministic behaviour 
can be found, recovering the rate equations of an equivalent PLM description in the 
K — > oo limit, and a Fokker-Planck equation (FPE) can be derived to characterise the 
fluctuations when K is finite. 

The importance of using the stochastic description of the dynamics when K is finite 
is revealed by the profound effect noise can have on the behaviour of the system. A 
striking example comes from the recent work of McKane and Newman p|] in the context 
of a predator-prey system. They show that although population cycles are absent in the 
deterministic limit of a simple system, there are sustained oscillations (cycles) in finite 
populations that can be explained by resonant amplification in the demographic noise. 
The origin of these cycles is a large-amplitude excitation of a deterministic spiral into 
the fixed point by the demographic noise in the system. In the language of Tome and 
de Oliveira [5], the demographic noise converts the deterministic, damped oscillations 
of the spiral, to undamped oscillations in the form of phase-forgetting quasicycles in the 
species densities. Subsequently, this phenomenon has been realised in models of systems 
as diverse as gene regulation [5] and measles epidemics [7]. 

These cycles also appear to be evident in (finite) Lotka-Volterra systems when 
the stochastic dynamics are embedded in real space (see e.g. [E]). The importance of 
this work is that it establishes a generic mechanism for sustained population cycles in 
predator-prey systems. The periodic orbits present in the deterministic Lotka-Volterra 
system, for example, are not robust to modifications of the dynamics or noise [2J: 
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intrinsic demographic noise has been shown to destabilise marginally stable predator 
prey cycles [9j [TUl [TT] . shifting the phase space trajectory between different limit cycle 
orbits until eventually one crosses an absorbing state where a species becomes extinct. 
The need for models in which cycles are robustly observed arises from the fact that 
some natural populations indeed exhibit cyclicity in abundances: for example, voles 
and lemmings exhibit complex multi-year cycles [12]. It is unlikely that such behaviour 
is achieved by fine-tuning of parameters with highly specific values that happen to favour 
cyclicity. 

In this work, we investigate a potential and hitherto unexplored mechanism for 
generating cycles in a population dynamical system without parameter tuning. It 
is based on the observation that when one has more than one stochastic variable 
(as is the case in a multi-species population dynamics) a nonequilibrium steady state 
that exhibits cyclic probability currents in configuration space typically arises. The 
question we address in this paper is whether these abstract cycles can be manifested as 
macroscopically observable cycles in species abundances. This is a reasonable hypothesis 
since it is known that in physical systems, such currents are directly observable in 
experiments on optically-tapped colloids [13J. There, the presence of a current was 
inferred from the persistent bias seen in the evolution of the polar coordinate of the 
colloid's motion. 

The fundamental origin of these probability currents can be understood from the 
Fokker-Planck equation [14] for the species frequencies Xj. This takes the form 

where the probability current J, can be expressed via the Kramers-Moyal expansion 
coefficients {a(x)} as 

d 

3 

One way to obtain a steady state is if currents vanish everywhere, J = 0. This 
corresponds to an equilibrium steady state in which the following conditions for detailed 
balance are met: 

• The drift term obeys the potential condition: dcti/dxj = dacj/xi ; 

• The diffusion term is thermal: a^- = const x 5ij. 

In a one-dimensional (ID) system — by which we refer to the size of the phase space; 
all the work presented here considers non-spatial systems — it is possible to make a 
change of variable such that these conditions are satisfied, unless a current is enforced 
on the system by virtue of periodic boundary conditions, or driving of the system at 
the boundaries. If the system has natural boundary conditions, where the probability 
and the current vanish at the boundaries, then necessarily the current must be zero 
everywhere in the steady state. In two or more dimensions, however, detailed balance 
is typically only satisfied if it is imposed from the outset, as is appropriate for systems 
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that are at thermal equilibrium. In particular, when one is dealing with systems that 
are defined by their dynamics, rather than by appealing to an energy function, a steady 
state in which the current is nonzero (but whose divergence vanishes) is most likely. 
In two dimensions (2D), the the vanishing of a divergence implies that 'streamlines' 
following the probability current form closed loops in the 2D phase space within the 
boundaries of the system, regardless of whether the system has natural, periodic or 
driven boundary conditions. 

The paper is organised as follows. In section [2] we examine the simplest 2D model 
of population dynamics in which a nonequilibrium current can be induced. This model 
has a stable fixed point, and the streamlines in probability space are described by 
ellipses. However, we show that these ellipses are not evident macroscopically. The 
simplest model that we have been able to find that exhibits cyclic behaviour purely 
as a consequence of nonequilibrium probability currents is set out in section [3j The 
deterministic limit of this model has dynamics that has a closed line (or manifold) of 
neutrally-stable fixed points. The dynamics are such that the system initially evolves 
towards the manifold. Once it is reached, one does not expect any further evolution in 
the deterministic limit: all forces vanish on the manifold. Naively, one would expect 
noise to generate unbiased diffusion around the manifold. However, the fact that detailed 
balance is not satisfied leads to a biased diffusion around the manifold. In section [5] we 
show that this current can be observed directly in simulations of the stochastic dynamics 
using the Gillespie algorithm. We furthermore show that a good approximation to these 
dynamics can be obtained using a reduction to a ID diffusion. Crucially, the noise in 
this approximation is multiplicative which, after transformation to a standard thermal 
diffusion problem, yields a complex non-conservative forcing. This dynamics can then be 
modelled as a biased diffusion between successive potential minima, the rates of which 
can be estimated using Kramers' escape rate theory. 

2. Probability currents near a stable fixed point 

The simplest population dynamical model with a nonequilibrium steady state has two 
species, A and B inhabiting a non-spatial patch of fixed finite size N = Xa + Xb + Xe- 
Here, Xa, Xb and Xe are the numbers of A individuals, B individuals and unoccupied 
sites. Individuals of both species reproduce with rate a, while each can die due to: (i) 
natural death with rate d; (ii) predation from the other species with rate p; and (iii) 
cannibalism from another of its species with rate c. Following the formalism for reaction 
kinetics of [15] , these processes give transition rates for the master equation of the form 

T(X A + 1, X B \X A , X b ) = 2a^(iV - X A - X B ) 
T(X A ,X B + l\X A ,X B ) = 2a^(N - X A - X B ) 
T{X A - 1, X B \X A , X B ) = dX A + c ^ - 1} + 2p^ 
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T(X Al X B - 1\X A ,X B ) = dX B + c Xb{X ° ~ 1} + 2p X j^ . (4) 

The factors of 2 arise from the combinatorics of the reaction kinetics; the rate constant 
d has been rescaled by iV and a, p,c by N — 1, bringing each rate into the van Kampen 
scaling form, where they are proportional to N [T6J, [3] . 

We carry out a systematic system size expansion of the master equation a la van 
Kampen [3] using the ansatz 

X A = Np A (t) + y/Nx A (t) , X B = N PB (t) + VNx B (t) (5) 

where p AyB are well defined densities and x A>B are now the random variables. To leading 
order in N, the deterministic behaviour predicts a commonly shared stable fixed point 

= 2a -d 

P *~ {2p-c){Aa + 2p + c) ' [) 

To the next order in N we find the FPE for the random variables x AjB which 
describe the fluctuations about this fixed point. Writing the vector x = (x A ,x B ), it is 

d -^^ = -^(TxP) + D-^- (7) 

dt dx^ }+ ^dxidxj [() 

where the drift matrix 7 and the diffusion matrix D 

( 2a + c 2(a + p) \ 

1 = - p *\2{a + p) 2a + c (8) 




D = 2ap,- ^ '-[ n I (9) 

^ Aa + 2p + c { M ' ' 

are independent of x. A generic feature of the van Kampen expansion is that this FPE 
has a drift term linear in x and a diffusion matrix that is constant. This implies that 
its steady-state solution always has the Gaussian form [3] 

P 5 (x) = - e -* XiS w (10) 

Zj 

where the inverse covariance matrix S satisfies the constraints 

7 T S + S 7 =2SDS (11) 

Tr(DS - 7) = . (12) 

The Gaussian form leads to the steady-state probability current vector J taking 
the form 

J 4 = (DS - l)i 3 x 3 Ps . (13) 

Expressing the probability current J = xPs where x is the average velocity vector and 
comparing to (13) we have for the fluctuations that 

x = (DS - 7)x = Ifx . (14) 

For the transitions rates given in Q we find that W = 0, meaning the system is in 
thermal equilibrium. From inspection of (JsJ) and (J9J) we see this is because the two 
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Figure 1. Time evolution of the total average angular displacement (</>(£)) measured by 
simulation of the stochastic dynamics Q by the Gillespie algorithm. The parameters 
are a — 2.1, d — 0.5, c = 2.5, p = 0.9. For the red line (top) r\ = —0.6, r 2 = 0, the 
black line (middle) r-y = r 2 = 0, the green line (bottom) r\ = 0, r 2 = —0.6. Each is 
averaged over 100 runs. 



necessary conditions for detailed balance are met: the potential condition is satisfied by 
7i2 = 721 and the diffusion is thermal, Dij = D eq 8ij . 

We can however induce a nonequilibrium current without modifying the 
deterministic contribution to the dynamics by introducing the parameters r^ and 
transforming the reaction rates for the birth and death processes for species A (and 
similarly for B) by 

a — >• a — 7"i , d — > d — 2?r , p ^- p + r 1 , c — > c + 2r\ . (15) 

Then, 7 remains constant, but the noise is modified so that 

D tj = [D eq + 2r iP *(2p* - 1)] 5^ (16) 

meaning that if r\ 7^ r 2 then D n 7^ D 2 2, making the noise athermal. This is sufficient 



to render the matrix W in (14) non-zero, thereby generating a current. The eigenvalues 



of W, which can be expressed as 



A. 



TrW v /(Tr^) 2 -4detVT 



(17) 



are purely imaginary due to the constraint (12). This means the solution to the 



autonomous system ( 14 ) are ellipses with a frequency of orbit 



UJ E 



Vdet W . 



We calculate the determinant of W by constructing the Gaussian matrix S in (10) by 
explicitly solving the linear FPE Q by the method of characteristics [I 



We now investigate how such a nonequilibrium current is manifested in the 
macroscopic population dynamics. Our primary tool is simulation of the dynamics 
using the Gillespie algorithm [T7j. To identify a sustained cyclic current in the system, 
we define for any given transition the angular displacement 8(f) about the fixed point p* 
If we sum up all the 8(p over multiple transitions, yielding </>(t), we should find that its 
ensemble average grows linearly in time, (<p{t)) ~ uiQt. 
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(a) ri = r 2 = 



(b) n =-0.6,r 2 = 



Figure 2. The phase space of the species A and -B. The system diffuses noisily about 
the fixed point, which for a — 2.1, d = 0.5, c = 2.5, p = 0.9 lies at Xa = Xb = 291. 



We see this is borne out in figure [TJ where the criteria for a current flowing are 
fulfilled through the parameters r^. For the parameters quoted in the figure we find 
by linear regression that the average angular velocity has a magnitude of loq = 0.97 for 
the (red) top and (green) bottom plots. We also see that the direction of circulation 
is determined by the athermality of the noise, i.e. if _D n is larger or smaller than D 22 . 
This measured speed tallies very well with the frequency of the elliptical orbit calculated 
by (18) for the same parameters: uje = 0.98. 



A crucial feature of the dynamics is that one does not observe cyclicity in population 
abundances in single realisations of the stochastic dynamics on short timescales. More 
precisely, it is not possible to see the precession of a single orbit over a 2tt interval. 
This is because rather than follow a fixed orbital path, the system diffuses noisily about 
the fixed point, meaning the cases of the system reaching thermal equilibrium or a 
nonequilibrium steady state are barely distinguishable, as is shown in figure [2j 

Therefore we find that a nonequilibrium current could not be advanced as an 
explanation for cyclic behaviour in natural populations, as observed in single field 
experiments. Moreover, our definition of an angle about the fixed point is tenuous 
as the system can diffuse very close to and even jump over this origin, making the 
notion of an orbit implausible. We are therefore drawn to the question: Is it possible 
to set up a population dynamical model so that cyclic behaviour is observed over in a 
single trajectory over short times? 



3. Macroscopic cycles on a neutrally-stable circular manifold 

The basic reason that observable cycles were not evident in the simple model of Section [2] 
is that the dynamics were attracted deterministically to a fixed point. In order for a 
macroscopic cycle to be evident, what is needed is for the dynamics to be repelled 
from some origin, so that an orbit may then be set up. In particular, if the dynamics 
are constructed such that the deterministic forcing vanishes along some closed line (or 
manifold) in phase space, and such that the dynamics are attracted onto the manifold 
from outside, then it could in principle be possible for nonequilibrium fluctuations to 
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generate a biased diffusion on the manifold. If so, in general we would expect to 
observe cyclical behaviour about the circular manifold. We freely admit that, from 
the perspective of modelling real biological populations, the dynamics that have such a 
property are somewhat contrived. However, from the more fundamental perspective of 
nonequilibrium statistical mechanics, we believe that this is the simplest such system 
that admits macroscopically observable currents arising from nonequilibrium probability 
currents. It is also a natural extension of quasi-neutral models of population dynamics 
from open to closed neutral manifolds. In [18], Parsons and Quince study a model of a 
non-spatial population with two species, each of which undergo the same birth and death 
processes. The model is non-neutral in the sense that the birth and death rate differ 
for each species. However, when the ratios of these two parameters are the same for 
both species, deterministically both population densities evolve to reside at a neutrally 
stable linear manifold, made up of a continuum of fixed points. In this sense the model 
is said to be quasi-neutral. 

This model is defined as follows.lt has a similar structure to that of the previous 
section, except that the size of the non-spatial patch is no longer conserved, and so the 
volume of empty space does not explicitly enter into the birth and death rates. Defining 
the intensive population numbers x A = Xa/K, xb = Xb/K where K is the carrying 
capacity of the system, both species undergo the following birth and death processes: 

T(X A + 1, X B \X A , X B ) = bx A (2 + 3x 2 A + x\ + 2x A x B ) birth of A 

T{X Al X B + \\X A , X B ) = bx B (2 + 3x 2 B + x\ + 2x A x B ) birth of B 

T(Xa — 1, X b \Xa, Xb) = bxA (q + (4 — q)xA + 2xb + x a + xax b ) death of A 

T(X A , X B - 1 \X A , X B ) = bx B (q + (4 - q)x B + 2x A + x B + x B x\) death of B (19) 

where b = b(K) and q are constants. As will be shown shortly, these rates have been 
designed to generate a circular manifold of neutrally-stable fixed points. To interpret 
them in terms of biological processes, we note that a term in x^Xg 1 corresponds to an 
interaction between n individuals of species A and m individuals of species B. So, for 
example, the birth of A at a rate proportional to x A would arise from some interaction 
between three individuals of species A, i.e., some cooperative interaction. Likewise, a 
birth rate for A proportional to xax b would arise through some interaction between 
individuals of species B: for example, the cooperative production of some resource by 
individuals of B that is beneficial to A. The higher-order terms in the death rates serve 
to stop the population sizes running out of control: they can therefore be interpreted 
as some kind of resource depletion implied by large populations. Thus although the 
combination of processes that yields a circular neutrally-stable manifold is somewhat 
specific, the biological principles involved (cooperative behaviour and resource depletion) 
are not entirely unreasonable. 

For analysis purposes, the useful quantities obtained from the rates are the moments 
of Sxa and Sxb- These are given by 

((6xA)\6x B y)) = J2 (SxA)\SxByT(X' A ,X' B \X A ,X B )r (20) 
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where X' A , B = Xa/b ± 1 and T(X' A , X' b \Xa, Xb)t is the probability that xa and/or xb 
changes by ±1/K in a time r. We write the moments in terms of the extensive variables 
b and K, and the intensive variables xa and xb via the functions Mij(xA,xs)'- 

b 



{(6x A y(6x B y)) 



Ki+3 



rMi a . 



(21) 



We find for the first and second moments: 

M 1)0 = xa (2 — q — (4 — q)x A — 2x^ + 3x A + x\ + 2xa%b ~ x a — x a% 2 b) 

M 0>1 = xb (2 — q — (4 — q)xs — 2xa + 3x| + x\ + Zxa^b — x B — xbx 2 a ) 

M 2 fl = xa (2 + q + (4 - q)xA + 2x B + 3x^ + x\ + 2xax b + x^ + x^l) 

M ,2 = ^b (2 + q + (4 - q)x B + 2x A + 3x| + x\ + 2xax b + x| + x B a;^) 

Mi,i = . (22) 

The deterministic behaviour can be found from the rate equations for (xa) and (xb)' 



d(x A ) 

dt 
d(x B ) 

dt 



lim 
lim 

MO 



(x A (t + T))-(x A (t)) 



(x B {t 



T 
T)) 



(Mt)) 



T 



M lfi 



M „ 



(23) 



can be computed using (22). In the limit K — > oo where we neglect fluctuations, 
(xA{t)) = XA{t) and we find after some algebra: 



x A = bx A (x A -!)[<?- (x A ~ I) 2 
x B = bx B {x B -!)[?- (^ - I) 2 
Defining the polar coordinates r and via 

x' A = xa — 1 = r cos(^) 



(x B 
(«b 



,r 



b 



x b 



rsin(0) 



we find that the fixed points of (24) are 

(x A = 0,x B = 0) 
(x A = l,x B = 0) 
(x A = 0,x B = 1) 
(x A = l,x B = 1) 
r = r = y[q . 



(24) 



(25) 

(26) 
(27) 
(28) 
(29) 
(30) 



The fixed point (26) is an absorbing state at the origin where both species are extinct 



and is repulsive. The fixed point (29) is also repulsive and sits at the centre of the polar 



coordinate system defined in (25). Physically this means that (K,K) is the centre of 
the circle in the extensive system. Two of the other fixed points are saddle points which 
sit on the A and B axis respectively, which are absorbing for the respective species. The 



fixed point (30) is a circular manifold of fixed points, located at fixed radius r from the 



centre given by (29) and is neutrally stable. 



To summarise, we find that at the level of the deterministic equations, the system 
evolves to a fixed point on the circular manifold. The fixed point that is reached would 
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Figure 3. Typical dynamical trajectory of the system defined by X^(t) and Xg(t) 



and the transition rates given in (19 1. These data were obtained by simulation using 
the Gillespie algorithm of the dynamics with parameters K = 5000, b = K 2 , 



0.5. 



be determined by the initial conditions. For finite K however, we expect the system 
to evolve to the neighbourhood of the manifold, and then diffuse about it. A typical 
snapshot of the evolution of the system is given in figure [3] this indeed shows that the 
occupied region of phase space is a circular annulus of large radius relative to its width. 



4. Reduction to a one-dimensional diffusion 

To analyse the model introduced in the previous section, we reduce the full dynamics 
in the two-dimensional space spanned by xa and xb to the diffusion of the angular 
coordinate in a polar coordinate system defined in (25) whose origin is the centre of 



the circle. A similar approach was used by Parker and Kamenev in studying the stability 
of a stochastic Lotka-Volterra model [TU]. There, as the angular component of the 
motion relaxes rapidly, it can be integrated out of the probability distribution, allowing 
the dynamics of the 2D system to be described by the ID stochastic radial motion 
between the deterministic limit cycles. In contrast, our main assumption is to neglect 
diffusion in the radial r direction off the stable circular manifold. That is, we assume 
that the restoring force that acts perpendicular to the manifold is sufficiently strong 
that any deviation away from the manifold does not contribute to the dynamics in an 
important way. We note, however, that this lateral diffusion has been seen to enter into 
an effective description on an open manifold [18] . Here, our aim is to see how well we can 
understand the full two-dimensional diffusion, and in particular any sustained angular 
velocity, within a highly simplified approximation. In a very recent piece of work [19] , 
Constable et al studied the dynamics of a two-species population which deterministically 
evolves to exist on a stable hyperbolic manifold. They carried out a similar dimensional 
reduction to that presented here and found it to be a good approximation of the full 
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(a) M 



(b) g{4>) 



Figure 4. The drift term / and the diffusion term g of the Langevin equation ( 34 ) 
for q — 0.5. 



2D dynamics. However a crucial difference is that the hyperbola has natural boundary 
conditions. As mentioned previously for a ID system, natural boundary conditions 
means the current must vanish at the boundaries. Therefore the system cannot support 
a sustained constant current in the steady state as it must be zero everywhere on the 
hyperbola. 

In order to study the diffusion in the <fi direction around the manifold, we derive a 
polar FPE for the evolution of the probability density P(<f),t), details of which can be 



found in Appendix A It is 

dP{<j>, t) __ d 

Ft " ~d6 



f(<t>)P 



+ 



1 d 2 



sPMP 



(31) 



where the drift term / and the diffusion term g are 

/(<£> = A ( 2r o( 12 + r o 2 )[sin(0) - cos(0)] + 4r (6 + r o 2 )[sin(30) + cos(30)] 
4r v 



+2Oosin(40) + 2rjj[sin(50) - cos(5</>)] 



(32) 



K 



64 + 28r 2 + r (56 + 6r 2 )[sin(0) + cos(0)] + 4r (6 + r 2 )[sin(30) - cos(3^)] 



(33) 



+24r 2 sin(20) - 20r 2 cos(40) - 4r 2 [sin(50) - cos(5^)] 
We can write the equivalent Langevin equation 

= /(0) + <?(0)^ (34) 

using the Ito prescription [H] where %)$ is Gaussian white noise with zero mean and unit 
variance. As one would expect, / and g are 27r-periodic as can be seen in figure Q. 



A key feature of the Langevin equation (34) is that the noise is multiplicative. The 



criteria for a nonequilibrium steady state to exist in a system with a FPE of the form 



(31) with periodic boundary conditions are more easily checked in the case where the 



noise is additive. It is always possible to perform a change of variable which transforms 
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a ID FPE with athermal diffusion to one with thermal diffusion [14] . The required 
change of variable is 



'2D 



o g(V) 

under which the FPE (31) becomes 



(35) 



d>=4>{9) 



(37) 



d t Q(6, t) = -d e (F(9)Q) + DdlQ (36) 

describing the time evolution of the probability density Q(6,t) where now the diffusion 
is a constant D and the drift force F(9) is 

The corresponding Langevin equation (a la Ito) is 

6 = F(9) + V2D r)e . (38) 

In this transformation one can choose D arbitrarily. For completeness we will proceed 
with arbitrary D, but in numerical calculations will set it to unity. 

Unfortunately, it is difficult to obtain an analytic expression for F(6), as one needs 



to evaluate the integral (35) and invert the resulting expression to obtain the function 



(j){9) that appears in (37). It is however possible to determine whether the drift force 
F{6) is conservative, and hence infer whether the system will reach thermal equilibrium, 
without an explicit expression for this function. 

First, we note that F(9) is conservative if the work done in one circulation of the 
manifold vanishes, i.e., if 

>d9F(9) = 0. (39) 

We can write an ansatz for F: 

F(9) =lu- d e V(6) , (40) 

where the first term wisa constant, non-equilibrium driving force and the second term 
is a conservative force derived from a potential V(8). To determine whether u — 0, 
formally, we should be able to write F{9) as a Fourier series 

'n9\ , (n6 s 

— I + b n sin I — 



1 oo 

no) = ^o + J2 



n=l 



a„ cos 



(41) 



where, using (35) 



'2D 



(42) 



2 J-«gtf) 

is the half-width of the transformed interval [— ir, n], viz, T = 8(tt) — 6(0). (We can take 
8(0) = with no loss of generality). The coefficients are given in the usual way as 

a n =- I d9F(9)cos' 

b n = — . .....,..,...., 

T / rr W V T 



T 



-T 

d6 F(6) sin ' 

T 



(43) 
(44) 
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Comparing (40) and (41) we see that u = a /2. So the question of whether F{6) is 



conservative and hence whether the system is reaches thermal equilibrium is equivalent 
to finding if cto is zero. 
Taking 

00 = ^1 dOF(0) (45) 

we apply the change of variable ( p5J ) : 

1 /" ?i(T) 6.6 



a 



(46) 



By definition 0(±T) = ±7r, F (6 [(/)]) can be obtained directly by (37), formally via 
(f){9 [(p\) = 0, and we know the Jacobian of the transformation. So, 



«o 



/ID 

/2D 



( fW d , u\ 

Uw-dV 11 ^ 



(47) 
(48) 



since the second term vanishes due to the periodicity of g(<fi). 

Although we have successfully side-stepped the problem of evaluating 0(6*), it 



remains the case that the integral (48) does not have a convenient closed form. However 



we can determine whether it is zero by finding if the integrand 



h(<f>) 



g 2 {4>) 



(49) 



is odd over the limits of integration. In figure [4] we see that / is odd and g is even about 
(f> = 7r/4. This implies that h is also odd about 7r/4, namely: 

"K 

.2 ~ 



h((f>) = -h 



(50) 



Due to the 27r-periodicity of h we can write (48) as 

a = d4> h((f>) = / d(f> h(4>) = / d(p h{4>) + 
Changing variable to /3 — n/2 — <p in the second integral we have 



4 



d(f> h((f>) . (51) 



«o 



d<p h{<f>) + 



4 d/3M|-/3) 



(52) 



Now applying (50) we see that a = 0. 



Intuitively it seems correct that the system reaches thermal equilibrium as the 



model defined by the processes in (19) is neutral and each species undergoes the same 



processes at equivalent rates. The insight that is gained from this analysis is that we now 
know in order to reach a nonequilibrium steady state, we must introduce processes which 
will enter the drift term / and the diffusion g in such a way as to make the integrand 
in (48) not be odd. As in Section [2j we are interested in the effect of different forms 
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of the noise whilst holding the deterministic contribution to the dynamics constant. To 
this end we introduce the following extra rates to our model: 

T(X A + l,X B \X A ,X B ) = b{p x x A x 3 B +p 2 x 3 a x B ) 
T{X A ,X B + l\X A ,X B ) = b{ Pl x A x 3 B +p 2 x 3 a x B ) 
T{X A - l,X B \X A ,X B ) = b(p 1 x A x 3 B +p 2 x 3 a x B ) 
T(X A ,X B - 1\X A ,X B ) = b(p 1 x A x B +p 2 x 3 a x B ) . (53) 

These new rates cancel in the first moments (x A ) and (x B ) and so leave the deterministic 



equations (23) unaltered. The second moments accrue the extra terms 



b 

((5x A ) 2 ) = —2 2 Pl x A x 3 B + 2p 2 x\x B r (54) 

((Sx B ) 2 ) = —2 2p 1 x A x 3 B + 2p 2 x A x B t . (55) 

iV o 



Substituting these into (32) leaves / as it was. However the form of g will change. 



The two new parameters p\ and p 2 will determine the magnitude and direction of the 



probability current. The new terms that appear in our expression for g (33) are 

2pi sin 2 (0)(l + r cos(0))(l + r sin(0)) 3 + 2p 2 cos 2 (0)(l + r sin(0))(l + r cos(0)) 3 . (56) 

From this we see that g will only remain even about 7r/4 as long as p\ = p 2 . Therefore 
as long as p\ ^ p 2 a probability current will flow in the steady state. This condition 
breaks the neutral selection of the model as the rates of birth and death for each species 
are no longer exactly equivalent. However quasi-neutrality is maintained in the sense 
that the deterministic behaviour of both species is the same [18] , each evolving to reside 
on the common circular manifold. 

5. Methods for measuring the current 

We now outline three approaches to measuring the current in the nonequilibrium steady 
state. The first is to use Monte Carlo simulation of the dynamics for the full problem. 
The second is to numerically integrate the Langevin equation within the one-dimensional 
reduction, thereby revealing any error that is introduced in this procedure. The third 
is to appeal to Kramers' escape-rate theory to estimate the current. Our expectation is 
these methods trade accuracy for precision and (in the latter case) analytical insight. 

5.1. Monte Carlo simulation 

The probability current can be written as 

J{<f>,t) = P{<f>,t)u (57) 

where P is the probability density and u is the average angular velocity. As already 
mentioned, we infer the flow of a current from a non-zero u. Using the Gillespie 



algorithm the stochastic model as defined by the rates in (19) and (53) can be simulated. 
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(a) Simulation 



(b) Numerical Integration 



Figure 5. Plot of the time evolution of the average aggregate angle obtained from: 
(a) simulation of the 2D stochastic model using the Gillespie algorithm; (b) numerical 
integration of the ID Langevin equation. For each, the parameters are K — 5000, 
b = K 2 , q = 0.5, p\ = 20 and -pi — 0, and the average is taken over 100 runs. 



We measure the total angular displacement </>g(£) which is positive in the anti-clockwise 
direction and quantifies the total distance travelled. For each update 

G (£ + r) = G + 50 (58) 



where 8(f) is calculated using (A.3). In figure 5(a) we plot the angular displacement and 
find that 

(4> G {t)) ~u G t. (59) 

For the parameters quoted, linear regression of the data yields uc = —0.97. 

5.2. Numerical integration of the Langevin equation 

The average velocity can also be calculated from direct numerical integration of the 



quasi-lD Langevin equation (34). We do so using the integration scheme [H] 

<f) L {t + dt) = <j> L {t) + f(4> L )dt + g((j)L)Vdt ^ . (60) 

For the same set of parameter values of previously, we find an ensemble-averaged angular 



displacement shown in figure 5(b) This time we find that lul = —1.27, in reasonable 



agreement with the simulations of the full 2D diffusion. (See below for a more detailed 
discussion of the different methods for estimating the current). 



5.3. Kramers' escape-rate theory 

Our final approach makes use of the transformation of the diffusion with multiplicative 
noise to a diffusion with additive noise described in Section HJ This allows us to estimate 
the nonequilibrium current via calculations of escape rates over potential barriers as done 
by Kramers |20j . following closely the presentation of the method in Hi 
The potential $(6*) that we are required to calculate is 



*(0) 



d6'F{6') 



(61) 
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(a) Pl = 20, P2 = 



(b) Pl = 0, p 2 = 



(c) P2 = 0, p 2 = 20 



Figure 6. $(#) for different values of the noncquilibrium parameters p\, p 2 with 
9 = 0.5. 



From earlier discussion, one would expect $(#) to only be a true potential if F{9) is 



conservative. However, using the ansatz (40) in (61) we have 

$(0) = -u e 6 + V(6) (62) 

up to an irrelevant constant. We know V(9) is periodic as it can be expressed as a 



Fourier series (41). Therefore in the domain of 9, V will be monotonically shifted by 



uo9, which we interpret as our potential $(#). 

Given the form of / and g derived for this model it is not possible to perform 
the transformation from (f) to 9 in (35), or compute the integral in (61) to obtain $. 



By approximating the function g with a mathematically tractable function we are able 
to make some progress in analytically approximating these required expressions. The 



technical details of this can be found in Appendix B 



In figure p\ are typical examples of &(9) according to different values of the non- 
equilibrium parameters p\ and p%. The thermal equilibrium condition p% = p2 is shown 



in figure 6(b) In this case there is no current flowing and the potential $(0) is periodic. 



However in the nonequilibrium cases shown in figures 6(a) and 6(c) we see that $(#) is 
not periodic. Though as the domain of of the potential is only 9 E [—9 , 9q], as defined 
in Appendix B it is only important that the upper and lower limits of the domain are 



the same to satisfy the periodic boundary conditions of the system. The absolute value 
of a potential is not important for deriving the force for it, so we do not require <&(9) to 
be periodic to regard it as being a potential. 

Each of the potentials in figure [6] are characterised by having a minimum between 
two barriers. We can derive the escape rate of a particle over these barriers; here we 



will do so for the right side barrier of the potential displayed in figure 6(c) We denote 



the position of the minimum in the well by 8w, the top of the barrier by 9b and the 



next minimum to the right of the barrier by A. For the FPE (36) in a nonequilibrium 
steady state the probability current 

J = -d e ^e{9)Q s {9) - Dd e Q s (9) (63) 

can be written as 

d 



J = -De-* B ' D 



89 



,$e/D 



Qt 



(64) 
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Integrating between 9w and A we have 

■A 



J [ dd e^ /D = D [e^ Bw)/D Q s {6 w ) - e^ {A)/D Q s (A)} 
Je w 



(65) 



If the barrier height A$ = <&(##) — <&(9w) is much greater than the diffusion D then 
the particle is far more likely to be found in the well about 9w This means we neglect 
the second term in the square brackets in (65), giving 

D Qs(O w ) ^ {dw)/D 



J 



^e/D 



(66) 



The current can be expressed as the probability p of being in the well at $w multiplied 
by the escape rate r from the well. Taking 9\ < 9w < #2 to define the domain of the 
well we write 

p = [ \e Q s (6) . (67) 



From (64) the stationary distribution is 

Q S (Q) = Ne -*W D - J e -W 



dO 



, e 



<S>{0')/D 



D 



(6* 



As we are interested in the stationary distribution in the well we introduce Qs(@w) by 
eliminating the normalisation constant iV: 

Je -m/D 



Qs(e) = Q 8 {0w)e-MO-*W D - 



D 



de > e *(?)/D _ / && e He>)/D 



Using this we can now write 



V = Qs(0w)e 



H8w)/D 



e-2 



d Q e -*(0w)/D 



(69) 



(70) 



0i 



where the contribution from the square bracket terms in (69) is negligible when 
considering the contribution from the domain of the well. Combining this with (66) 
we can express the escape rate as 

1 n 1 



62 &e e -*w D [ A de e*W D . 



- =l = - I d9 e -*WA" / d9 e*W' u . (71) 

r J D Je 1 Je w 

We Taylor expand each integrand in the above expression, the first about 9w, the 

second about 9 b'- 

, , &'(0w) 



*(0) 



$( 



$( 



7 w 



7 B 



2 

W'(e B ) 



e^ 2 



7 W , 



0,* 



7 B 



(72) 
(73) 



Using these second order expansions, we can extend the boundaries of each integral to 
±00, meaning that to both integrals in (71) are now Gaussian. Computing them leads 

<m,)i<i>"(M <-**/ D . (74) 



to the final expression for the escape rate 
1 



2n 
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Figure 7. Overlay of the different measurements of lu using the three distinct methods 
of section [5] plotted against a relative measure of the nonequilibrium current pi — Pi- 
For each data point either pi or p 2 is zero. 



To find the value of these escape rates, we numerically compute the value of the second 
derivatives at the minima and maxima after explicitly constructing the potential $(#) 
using the mathematical software Maple (TM). 

Denoting by r~ and r + the escape rate for the left and right barrier respectively, 
we express the average change in 9 due to hopping over the right or left barrier in a 
time t as 

(60) = (r+A e - r-A e )r (75) 

where A# = 2# is the distance between wells, i.e. the period of the system. In the limit 
t — > this gives the average angular velocity 

u K = A e (r + - r~) . (76) 



It is clear from ( 74 ) that a current exists due to the difference in the height of the two 



barriers. When the system is in thermal equilibrium as in figure 6(b) there is no current 
as we are equally likely to hop left or right around one circuit. With the parameter set 
used before, and setting D — 1, we find u>k = —0.83. Again, this value is in reasonable 
agreement with those previously obtained by other means. 



5.4- Comparison of the three methods 

To better understand how well the currents obtained via these different approaches 
correspond with each other, we compare in Figure [7] these measures for different values 
of the p parameter which controls the extent to which detailed balance is violated. We 
see that all three measurements obey the same qualitative trend, and remain within an 
order of magnitude of each other. 

This validates our principal approximation to reduce the full two dimensional 
system to one dimension, that of the polar angle 0, by neglecting radial diffusion. In 
particular, we note that the ID criterion for a non-equilibrium steady state pi = p 2 
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i 




Figure 8. A single realisation of 4>c obtained by simulation of the stochastic dynamics 
with parameters: K = 5000, b = K 2 , q = 0.5, p\ = 20, p 2 = 0. 



is borne out by simulations of the full 2D dynamics. The difference in the dynamic 
measurements, ujl is larger than uq, is understandable as while we neglect any radial 
diffusion in the ID treatment, it is still present in the simulations, as witnessed in 
figure [3j We expect time spent diffusing radially to slow the rate of polar diffusion. 

The main assumption of the Kramers escape rate calculation is that the ratio of 
the barrier size A$ to the diffusion constant D is very large. In practice, we set D — 1, 
and it is not possible to tune the model to allow us to independently control the barrier 
height and the diffusion strength. Typically, we find that the ratio is 2 < AQ/D < 3, 
which is a likely cause of the quantitative discrepancy between this method of obtaining 
uj with the other two. 

6. Discussion and Conclusions 



In this work, it has been our aim to understand the conditions under which a 
nonequilibrium steady-state current in a multi-species population dynamics model may 
be visible macroscopically in the form of population abundance cycles. The simplest 
candidate for such a system has an attractive fixed point and an absence of detailed 
balance. Here we found that cyclicity is not evident on the timescale of a single 
(probability) cycle in a single realisation of the dynamics. One must either perform an 
average over long times or an large ensemble to uncover this weak systematic cyclicity. 

One situation that allows these currents to be observed in single realisations of 
the dynamics has a deterministic dynamics that is neutral on some closed manifold 
in configuration space. By this we mean that all deterministic forces vanish on the 
manifold, and that the manifold itself is attractive. A nonequilibrium steady state can 
then be manifested as a systematic orbit around the manifold. Although, the specific 
interactions that yield such a structure are somewhat contrived from a purely biological 
viewpoint, this work highlights some interesting features of stochastic dynamical systems 
from the perspective of nonequilibrium statistical mechanics. 

Most notably, the stochastic trajectories along the manifold are somewhat complex, 
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as figure [8] shows. We see that the system tends to diffuses over a small region of the 
phase space before sharply jumping a distance of about w 2tt radians, i.e. a full circuit 
of the manifold. The origin of this motion can be understood from the potential picture 
after mapping to additive (thermal) noise as described in section |4J After transforming 
the multiplicative noise, the deterministic equations acquire additional terms that can 
be interpreted as a constant driving force acting in a periodic potential with multiple 
maxima and minima on the manifold. The dynamics will reside for some time in a 
potential minimum before escaping over one of the barriers to a neighbouring minimum. 
Since the potential is periodic, one eventually returns to the same minimum (hence the 
2-7T jump). The driving force leads to an asymmetry in the barrier heights, which in turn 
yields a systematic current in one direction around the circle (i.e., a 'cycle', albeit not 
a strictly periodic one). 

This analysis leads us to believe that in general the closed manifold does not have 
to be neutrally stable in order to observe this cyclical behaviour, as long as the manifold 
itself is attractive from the outside. In such an instance, the non-zero angular velocity 
generated by the probability current would still drive the system over potential barriers, 
but now they barriers are due to the deterministic forces acting on the manifold. While 
the time scale for passing over a barrier would increase markedly, quasicycles would still 
be generically observable in a system with athermal noise whose deterministic dynamics 
evolve to a closed manifold. 

The observation of quasicycles due to a noise- induced non- conservative driving force 
highlights the importance of the form of the noise that is manifest in these stochastic 



population models. In the Langevin equation (34) the non-equilibrium processes given 



in (53) only appear in the diffusion term g. Therefore if one naively assumed the noise 



to be thermal, i.e. g is constant, integrating up the drift force / would yield a periodic 



potential similar to the one displayed in figure 6(b) , and one would conclude no current 
is flowing. We see that stochastic effects alone are responsible for a current to flow in 
the system, keeping the system out of equilibrium. 

In the context of modelling population dynamical systems more generally, our 
findings further showcase the advantages of taking an IBM approach to modelling 
finite population dynamics. The deterministic and stochastic contributions to the 
dynamics can be derived from the defined stochastic processes of the model, so one can 
analytically understand and also control the effects of noise. We have found that finite 
size populations in which the dynamics are non-neutral, never relax fully to equilibrium, 
but instead inhabit a steady state where their is a thermal bias in the fluctuations due 
to the presence of a non-equilibrium probability current. 

In this work we resorted to a number of approximations to compute the steady state 
current in the dynamical system. The first of these was a reduction to a single-coordinate 
description by disregarding one of the degrees of freedom in the system. Sophisticated 
methods have been applied to integrate out this degree of freedom in the context of quasi- 
neutral diffusion along an open interval [T5] (as opposed to one that is closed/periodic, 
as here). It would be of interest to see if similar methods can be applied to determine 
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what is lost in such a dimensional reduction, as this may be of utility in understanding 
high-dimensional stochastic dynamical systems more generally. Moreover, we made 
various approximations in order to apply Kramers' escape rate theory to diffusion on 
the manifold: it may be that more direct approaches to estimating the current in such 
systems can be found. Finally, and more generally, it would be interesting to establish 
if there are other ways that a nonequilibrium current may enter into the macroscopic 
dynamics of a stochastic population dynamical system in ways that are not immediately 
evident from their deterministic counterparts. 
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Appendix A. The polar Fokker-Planck equation 



To derive a Fokker-Planck equation for the evolution of the probability density P(<p, t) 
we terminate the Kramers Moyal forward expansion [1] 



dt 



n=l 



77! 



a n {<j>)P{(j>,t) 



(A.l) 



at the second term. This truncation is valid if ol\ and a 2 are the only non-zero jump 
moments as defined by 



a, 



lim 

T->0 



([0(t + r)-0(t)H 



lim 

T->0 



((<¥)"} 



(A.2) 



T r->0 T 

The moments of 5<f) can be obtained from the moments of 6~xa and Sxb by the relation 

5(f)(6x'5x'\x'x') = tan" 1 ( X 'f + 5 fA - tan" 1 ( ^] . (A.3) 

\x A + 5x A J \x A J 

Taylor expanding to second order in 8x' A and Sx' B about 5x' A = 5x' B = and keeping 
terms up to order 5 2 yields 

X B ($: I \ i x a , 



- -3^'a) + dti^E) 



+ 



x A x B 



{8{x' A f - {6x' B f) + 



lx 



i \2 
BJ 



X 



I \2 



Applying the polar transformation (l25|) and averaging yields 



-(Sx' A Sx'2 



BJ 






b 

b 



sin(20) 



M 2>0 - M 0j2 



sin 2 (0)M 2>o + cos 2 (0)M o ■ 



T . 



(A.4) 

(A.5) 
(A.6) 



Probability currents in finite populations 



22 



where we use the fact that deterministically we are at a fixed point so (5xa) = (8xb) = 0. 



For the higher order jump moments we can infer from (21 ) and rt22|) that to leading 
order in K 



({8<f>) n )~bK- n G n {x A ,x B )T, n>3 



(A.7) 



where G n are functions of the intensive variables, independent of K. Using (A.2), we 
have from the above expressions for the moments that 

b 



cti,ai2 



O; 



lim 

&->oo K 2 



lim 



K' 



n > 3 



(A.8) 
(A.9) 



The change in limit is an observation that the rates defined in ( 19 ) scale with b and so 



the time r until an event happens scales like 1/6. This means that in order to have a 
non-zero first and second jump moment with higher order jump moments vanishing we 
should choose b = K 2 . 



Appendix B. Approximating 9{4>) and <&(9) 
We are faced with the task of computing 



*(9) 



d9 ; F{9') 



-V2D 






:%(0) 



9') 



Applying the same change of variable as in (46) 

rflfl) 

$(0) = -V2D 



471) we have 



( f{4>) 

W(<t>) 



g{4>) 



(b.i; 



(B.2) 



(B-3) 



The form of g((f>) as it is given in (33 ) makes the above integral intractable. We therefore 
make the approximation G(<p) ~ G g ((p) = l/g(<f)) where G has the form 

G(0)=c G + A G cos(0-6 G ). (B.4) 

The phase 6 G is set to by matching numerically the position of the first extremum of 
G g and G for > 0. The other two parameters are set by making the maximum and 
minimum values of G„ and G the same. We illustrate this approximation for a typical g 



in figure Bl Applying the approximation G{<p) 

-V2D 



§o(9) 



d0 /(0)G(0) 



1/(7(0) our expression (B.3) becomes 

(B.5) 



l fyG(0) \ 
2 G(0) J ■ 

It is possible to numerically integrate this integral by computer 
evaluate it in the 9 coordinate we must have an analytical form of 
into it. 



However, to 
to substitute 
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Figure Bl. The inverted diffusion term 1/g (red, solid) and its approximation G 
(green, dashed). 



The transformation from 9 to in (35) does not have a closed form for the exact 



g{4>) derived for this model so we approximate it by 



0(0) « M0) 



4> 
''2D I d(f)'G{(t)') 
'0 



(B.6) 



which using (B.4) gives 



9 G {4>) = V2D(c G( p - A G sin(0 - b G ) - A G sin{b G )) . (B.7) 

This functional relation is shown in figure E |2(a)[ What we really require is the inverted 
form of this relation, giving us <fi in terms of 9. The inversion is not possible, however 
we can approximate it via 

d0 



9, 



'2D 



lQ c t + A t cos((j)) ' 
This integral has the solution 



9, 



2V2D 
Cty/1 - a t 



tan 




tan 




(B.8) 



(B.9) 



where q = At/ct- We need to approximate the values q and a t to give a good fit with 
the form of 9(<p) from (B.7). 

We do this by matching the linear gradient and the curvature near <fi = of 9 G (<p) 
and 9t(4>). The slope we match by finding the straight line gradient between two points 



that lie on 9 t ((j)) and making it equal to slope of the constant term in (B.7). Choosing 
the two points <fi = —n and <fi = tc we have from (B.9) that 

\[2D n 



9t(±n) 



Matching the gradients: 

6 t (n) - 6 t (-n) 
2% 



'2Dc G 



(B.10) 



(B.ll) 
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(a) 0gW 




(b) e G (<l>) and e t (4>) 



Figure B2. Comparison of the two approximations of the transformation 6g(4>) ( re d, 
solid) and t (4>) (green, dashed). 



yields 



1 



CG 



c t Jl -a? 



To match the curvature near the origin we rearrange (]B.9|) to 

Cty/1 ~ Of 9 t 



tan 



a, 



'2D 



>/l" 



tan 



Substituting 6q for t using (B.7) gives 
qa/1 -a? 6> G 



tan 



tan 



'2D 



Q 



vt 



2V2D 



/^(cg^ - A G sin(0 - bo) - A G sin(6 G )) 



Now setting the above two expressions equal and Taylor expanding about 
have 

l-O t <f) Cty/l -a\ ( .\ 

c G - A g cos(6g) 



T - a? 2 



which upon matching coefficients gives 

c G - A G cos(6 G ) • 



c t (l + a t ) 



(B.12) 



(B.13) 



(B.14) 
= we 

(B.15) 
(B.16) 



We can solve (B.12) and (B.16) simultaneously to find q and at for the approximate 



transformation (B.9). A comparison of the two approximate transformations are given 
in figure E |2(b)[ 



We can invert (B.9) and find the functional form we require: 

-l (cty/l-a%6 



l_-a| 
- Of 



tan 



'2D 



(B.17) 
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Finally, we need to find the limits of the domain in the 9 variable. Substituting the <fi 



limits ±7r into (B.9) we find 9 E [—9q,9q], with 



9 Q = ffL= . (B.18) 



Cfs/l-aj 
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